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I. INTRODUCTION AND SUMMARY 



It was shown that an all-sky full frequency bandwidth coherent search of data of many months duration for contin- 
uous gravitational-wave signals is computationally too prohibitive to be realized with presently available computing 
power • Consequently two alternative approaches emerged. One that we shall call tracking involves tracking of 
lines in the time-frequency plane built from the Fourier transforms of short (around 40 minutes long) stretches of 
data The other that we shall call stacking involves dividing the data into shorter (around a day long) stretches, 

searching each stretch for signals, and enhancing the detectability by incoherently summing the Fourier transforms 
of data stretches ||. In this paper we develop data analysis tools for the first stage of the stacking approach, namely 
for the coherent search of a few days long stretches of data. These tools involve storage of the data in the Fourier 
domain database, calculation of the precise position of the detector with respect to the solar system barycenter using 
JPL ephemeris, calculation of the detection thresholds, approximation of the signal by a simple linear model, and 
construction of the grid of templates in the parameter space ensuring no loss of signals. 

In this theoretical work we have in mind a particular set of data, namely the data collected by the EXPLORER bar 
detector. This detector is most sensitive over certain two narrow bandwidths of about 1 Hz wide at frequency around 
1 kHz. To make the search computationally feasible we propose an all-sky search of a few days long stretches of data 
in the narrow band where the detector has the best sensitivity. We find that with a 250 Mflops of computing power 
we can carry out a complete all-sky search in around a month for signals within a bandwidth of 0.76 Hz at frequency 
of 922 Hz and for observation time of 2 days. We assume that the minimum characteristic time for the change of the 
signal's frequency is 1000 years. By our search we shall be able to detect all the continuous gravitational-wave signals 
with the parameters given above and with the dimensionless amplitude larger than 2.8 x 10~ 23 with 99% confidence. 

The tools developed in this work can be applied not only to the analysis of the bar data but also to the interferometer 
data. In the case of the wide-band interferometer one can divide the data into many narrow bands and analyze each 
band using the algorithms presented in this paper. 

The plan of the paper is as follows. In Sec. II we present responses of both a laser interferometer and a resonant 
bar to a continuous gravitational wave. The optimal data processing method is derived in Sec. [II. In Sec. [V we 
describe the construction of the frequency domain database and in Sec. |v| we discuss how to calculate the precise 
position of the detector with respect to the solar system barycenter. In Sec. [y| we introd uce an approximate model 
of the detector's response called the linear model and we describe its performance. In Sec. VH we construct a grid of 
templates in the parameter space ensuring no loss of signals, we calculate the computational requirements to do the 
sear ch, an d we perform Monte Carlo simulations to test the performance of our grid and our search algorithm. In 
Sec. VIII we discuss and choose parameters for our planned all-sky search of the EXPLORER detector data. Several 
details are presented in the Appendixes. 
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Many methods presented in this paper rely on general analytic tools developed in a previous paper of this series || 
that we shall call Paper III. 

II. DETECTOR'S NOISE-FREE RESPONSE 

In this section we present the response of a laser interferometer and a resonant bar detector to a continuous 
gravitational wave. We include the effects of both amplitude and frequency modulation of the response. 

Dimensionless noise-free response function h of the gravitational- wave detector to a weak plane gravitational wave 
in the long wavelength approximation [i.e. when the size of the detector is much smaller than the reduced wavelength 
X/(2tt) of the wave] can be written as a linear combination of the wave polarization functions h + and h x ■ 

h(t) = F+(t) h+(t) + F x (t) h x (t), (2.1) 

where F + and F x are called the beam-pattern functions. 

A. Interferometric beam-pattern functions 

Noise-free response function h of the laser interferometric detector is defined as the diffe rence between the wave 
induced relative length changes of the two interferometer arms. Derivation of formula fl2.l| ) in the case of the laser 
interferometer can be found e.g. in Sec. II A of Ref. [Q. 

Explicit formulas for the interferometric beam-pattern functions F + and F x from Eq. ( |2.1[ ) are derived e.g. in 
Sec. II A of Ref. The functions read 

F+ (t) = sin C [a(t) cos 2ip + b(t) sin 2ip] , (2.2a) 
F x (t) = sin C [6(f) cos 2ip - a(t) sin 2ip] , (2.2b) 

where 

a(t) = - sin 27(1 + sin 2 </>)(l + sin 2 8) cos[2(a- <j> r - fi r f)] 

— — cos27sin0(l + sin 2 8) sin[2(a — <j) r — Q, r t)\ 
+ — sin 27 sin 2<p sin 28 cos(a — <j) r — Q r t) 

— — cos 27 cos cf> sin 28 sin(a — <p r — Q r t) 
3 

+ - sin 27 cos 2 0cos 2 8, (2.3a) 
6(f) = cos 27 sin (j) sin 8 cos [2 (a — (p r — fi r f)] 

+ - sin 27(1 + sin 2 </>) sin 8 sin[2(a — <fi r — Sl r t)] 
+ cos 27 cos <f> cos 8 cos(a — <p r — Q r t) 

+— sin 27 sin 20 cos 8 sin(a — 4> r — Q r t). (2.3b) 

In Eqs. ( |2.2| ) C is the angle between the interferometer arms (usually C = 90°) and ip is the polarization angle of the 
wave. In Eqs. (|2.3|) the angles a and 8 are respectively right ascension and declination of the gravitational- wave source. 



The geodetic latitude of the detector's site is denoted by 4> [hs precise definition is given in Eq. (5.1a) below], the 



angle 7 determines the orientation of the detector's arms with respect to local geographical directions: 7 is measured 
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counter-clockwise from East to the bisector of the interferometer arms. The rotational angular velocity of the Earth 
is denoted by f2 r , and r is a deterministic phase which defines the position of the Earth in its diurnal motion at t = 
(the sum r + Q r t essentially coincides with the local sidereal time of the detector's site, i.e. with the angle between 
the local meridian and the vernal point; see Sec. VII B below). 



B. Bar beam-pattern functions 



Derivation of explicit formulas for the bar beam-pattern functions F + and F x from Eq. (2.1) is relegated to Ap- 
pendix A of the present paper. The functions read 



where 



F + (t) = a(t) cos 2-0 + 6(f) sin 2^, 
F x (t) = b(t) cos 2-0 - a(t) sin 2-0, 

a(t) = -(cos 2 7 — sin 2 7sin 2 0)(1 + sin 2 5) cos[2(a — <p r — Q r t)] 



(2.4a) 
(2.4b) 



+— sin 27 sin 0(1 + sin 2 8) sin[2(a — <p r — fi r t)] 



- — sin 2 7 sin 20 sin 28 cos [a 



n r t] 



-— sin 27 cos sin 25 sin(a — r — fl r t) 



+ -(l-3sin 2 7cos 2 0)cos 2 (5, 

b(t) = — sin 27 sin sin S cos[2(a — r — f2 r t)] 

+ (cos 2 7 — sin 2 7 sin 2 0) sin(5sin[2(a — r — Cl r t)] 

— sin 27 cos cos 8 cos(a — r — fl r t) 

— sin 2 7sin20cos<5sin(a! — r — fl r t). 



(2.5a) 



(2.5b) 



In Eqs. (2.4) and (2.5) a, 8, ip, 0, f2 r , and r all have the same meaning as in Eqs. (2.2)— (2.3); the angle 7 determines 
the orientation of the bar detector with respect to local geographical directions: 7 is measured counter-clockwise from 
East to the bar's axis of symmetry. 



C. Wave polarization functions 



We are interested in a continuous gravitational wave, which is described by the wave polarization functions of the 
form 



h+(t) = h + cos*(t), 
h x (t) = h ax sin*(t), 



(2.6a) 
(2.6b) 



where h Q+ and ho x are the independent amplitudes of the two wave polarizations. These amplitudes depend on the 
physical mechanisms generating the gravitational wave. In the case of a wave originating from a spinning neutron 
star these amplitudes can be estimated by 



where / is the neutron star moment of inertia with respect to the rotation axis, e is the poloidal ellipticity of the star, 
and r is the distance to the star. The value of 10 -5 of the parameter e in the above estimate should be treated as an 
upper bound. In reality it may be several orders of magnitude less. 
The phase "J of the wave is given by 



(2.8a) 



si 
fc=0 



ffc+1 



n • r SSB {t) 



(fc+1)! 



^ k kV 



k=0 



(2.8b) 



In Eqs. (2.8) $o is the initial phase of the wave form, rgse is the vector joining the solar system barycenter (SSB) with 
the detector, no is the constant unit vector in the direction from the SSB to the neutron star. We assume that the 
gravitational wave form is almost monochromatic around some angular frequency luq which we define as instantaneous 
angular frequency evaluated at the SSB at t = 0, uik (fc = 1,2, . . .) is the kth time derivative of the instantaneous 
angular frequency at the SSB evaluated at t = 0. To obtain formulas (2.8) we model the frequency of the signal in 
the rest frame of the neutron star by a Taylor series. For the detailed derivation of the phase model (2.8) see Sec. 
II B and Appendix A of Ref. [|. The inte gers si and s 2 are the number of spin downs to be included in the two 
components of the phase. We need to include enough spin downs so that we have a sufficiently accurate model of the 
signal to extract it from the noise. This will depend on a given data and in particular on the length of the observation 
time. 



D. A linear representation 



It is con ve nien t to write th e resp onse of gravitational- wave detectors described above in the following form [cf. Eqs. 
O, (]2j), @, (||), and Q] 



i=l 

where the four constant amplitudes At are given by 

A\ = sinf (ho+ cos 2ip cos $o — ft-ox sin 2ip sin $o) > 

A 2 = sin£(/!.ox co s 2tp sin $ + h + sin 2-0 cos $o) > 
A 3 = sin £ (— ho+ cos 2\p sin 3> — ho x sin 2ip cos $0)7 

A4 = sin£ (ho+ cos 2ip cos $0 — ft-ox sin 2^ sin $0) ■ 

The amplitudes A4 depend on the physical mechanisms generating the gravitational wave. 
The time dependent functions hi have the form 



hi(t) 


= o(t) 


cos$(t), 


h 2 (t) 


= b(t) 


cos$(t), 


h 3 (t) 


= o(t) 


sin$(i), 


h 4 {t) 


= b(t) 


sin$(i), 



(2.9) 

(2.10a) 
(2.10b) 
(2.10c) 
(2.10d) 

(2.11a) 
(2.11b) 
(2.11c) 
(2. lid) 



where the functions a an d b are given by Eqs. (2.3) for an interferometer and by Eqs. (2.5) for a bar, and $ is the 
phase given by Eq. (2.8b). The modulation amplitudes a and b depend on the right ascension a and the declination 
S of the source (they also depend on the angles 4> and 7). The phase $ depends on the frequency uj , s spin-down 
parameters u>k (k = 1, . . . , s), and on the angles a, 6. We call parameters luq, w^, a, S the phase parameters. Moreover 
the phase $ depends on the latitude <p °f t ne detector. The whole signal h depends on s + 5 unknown parameters: 
/io+, h() X , q, (5, wo, Wfc. The formulas above apply to both resonant bar and laser interferometric detectors, but for the 
case of bars the angle £ is always 90°. The main feature of the above representation is that the signal is decomposed 
into a linear combination of several components. 
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III. OPTIMAL DATA ANALYSIS METHOD 



The signal given by Eq. ( |2.9| ) will be buried in the noise of a detector. Thus we are faced with the problem 
of detecting the signal and estimating its parameters. A standard method is the method of maximum likelihood 
(ML) detection that consists of maximizing the likelihood function, which we shall denote by A, with respect to the 
parameters of the signal. If the maximum of A exceeds a certain threshold calculated from the false alarm probability 
that we can afford, we say that the signal is detected. The values of the parameters that maximize A are said to be 
the maximum likelihood estimators of the parameters of the signal. The magnitude of the maximum of A determines 
the probability of detection of the signal. 

We assume that the noise n in the detector is an additive, stationary, Gaussian, and zero-mean continuous random 
process. Then the data x (if the signal h is present) can be written as 

x(t)=n(t)+h(t). (3.1) 

The log likelihood function has the form 



logA=(x\h)-±(h\h), (3.2) 



where the scalar product ( • | • ) is defined by 



o 



In Eq. (^^) tilde denotes the Fourier transform, asterisk means complex conjugation, and Sh is the one-sided spectral 
density of the detector's noise. 

Since we assume that over the bandwidth of the signal h the spectral density £/,(/) is nearly constant and equal to 



Sh(fa), where /o is the frequency of the signal measured at the SSB at t — 0, the scalar products in Eq. (3.2) can be 
approximated by 

(x\h)&-^— [ ' x(t)h(t)dt, (3.4a) 

Jh(Jo) Jo 

(h\h) « ^j- f° [h(tf dt, (3.4b) 

bh{jO) Jo 

where T is the observation time, and the observation interval is [0, T \. It is useful to introduce the following notation 

1 " T ° 



— \ x(t)dt. (3.5) 

J o JO 



Applying this notation and making use of Eqs. (3.4), the log likelihood ratio from Eq. (3.2) can be written as 



1 A 2T ' 

In A 



t) { {XH) ~ ^ ■ ^ 



Shift 

In Sec. Ill of Paper III we have analyzed in detail the likelihood ratio for the general case of a signal consisting 
of N narrow-band components. Here we only summarize the results of Paper III and adapt them to the case of our 



signal (2i)). The signal h depends linearly on four amplitudes Ai. The likelihood equations for the ML estimators A4 

0, * = !,... ,4. (3.7) 



of the amplitudes A^ are given by 

<91nA 



dAi 

One can easily find the explicit analytic solution to Eqs. ( |3.7| ). To simplify formulas we assume that the observation 
time T is an integer multiple of one sidereal day, i.e., T a — n(2Tr/Q r ) for some positive integer n. Then the time 
average of the product of the functions a and b vanishes, (ab) — 0, and the analytic formulas for the ML estimators 
of the amplitudes are given by 
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li«2^, (3.8a) 
(a 2 ) 

£«2<^>, (3.8b) 

1 3 « 2<^> , (3.8c) 

£«2<|^>. (3.8d) 

Explicit formulas for the time averages (a 2 ) and (b 2 ) are given in Appendix B. The reduced log likelihood function 
T is the log likeliho od f unction where amplitude parameters Ai were replaced by their estimators Ai. By virtue of 



Eqs. (3.8) from Eq. (3.6) one gets 



where 

F a := [ x(t)a(t)exp[-i$(t)]dt, (3.10a) 
•A) 

sc(i) b(t) exp[-i$(t)] dt. (3.10b) 

We thus see that the ML detection in Gaussian noise leads to a matched filter. The ML estimators of the signal's 
parameters are obtained in two steps. Firstly, the estimators of the frequency, the spin-down parameters, and the 
angles a and S are obtained by maximizing the functional T wit h re spect to these parameters. Secondly, the estimators 



of the amplitudes Ai are calculated from the analytic formulas (3J3) with the correlations (xhi) evaluated for the values 
of the parameters obtained in the first step. Thus filtering for the narrow-band gravitational-wave signal requires two 
linear complex filters. The amplitudes Ai of the signal depend on the physical mechanisms generating gravitational 
waves. If we know these mechanisms and consequently we know the dependence of Ai on a number of parameters we 
can estimate these parameters from the estimators of the amplitudes by least-squares method. 

We shall now summarize the statistical properties of the functional T. We first assume that the phase parameters 
are known and consequently that T is only a function of random data x (the case when the phase parameters are 
unknown and need to be estimated will be considered later in this section). When the signal is absent 2T has a \ 2 
distribution with four degrees of freedom and when the signal is present it has a noncentral x 2 distribution with four 
degrees of freedom and noncentrality parameter equal to the optimal signal-to-noise ratio d defined as 



d:= y/(h\h). (3.11) 

This is the maximum signal-to-noise ratio that can be achieved for a signal in additive noise with a linear filter [ . 
This fact does not depend on the statistics of the noise. Consequently the probability density functions po and p\ 
when respectively the signal is absent and present are given by 

Po(T) = Texp(-T), (3.12a) 

Pl (d, T) = ^-h {dV0) cxp (-? - ^d 2 ^j , (3.12b) 

where Ii is the modified Bessel function of the first kind and order 1. The false alarm probability Pp is the probability 
that T exceeds a certain threshold T Q when there is no signal. In our case we have 

/■oc 

P F {T ) := Po (T) dT=(l + T ) exp(-T ). (3.13) 



G 



The probability of detection Pjj is the probability that T exceeds the threshold T when the signal-to-noise ratio is 
equal to d: 

poo 

P D {d,T ):= J Pl {d,F)dF. (3.14) 

The integral in the above formula cannot be evaluated in terms of known special functions. We see that when the noise 
in the detector is Gaussian and the phase parameters are known the probability of detection of the signal depends on 
a single quantity: the optimal signal-to-noise ratio d. 

When the phase parameters are unknown the optimal statistics T depends not only on the random data x but also 
on the phase parameters that we shall denote by Such an object is called in the theory of stochastic processes a 
random field. Let us consider the correlation function of the random field 

C&O := E {[Hi) rn{Z)][F{ti') Mt)}} , (3-15) 

where 



(3.16) 



When the correlation function C depends only on the difference £ — the random field T is called homogeneous and 
for such fields we have developed in Paper III a theory to calculate the false alarm probabilities. 

The main idea is to divide the space of the phase parameters £ into elementary cells which boundary is determined 
by the characteristic correlation hypersurface of the random field T . The correlation hypersurface is defined by the 
requirement that at this hypersurface the correlation C equals half of the maximum value of C. Assuming that C 
attains its maximum value when £ — £ =0 the equation of the characteristic correlation hypersurface reads 



C(t) = \C(0), 



(3.17) 



where we have introduced r := £ — £ . Let us expand the left-hand side of Eq. (3.17) around r = up to terms of 
second order in r. We arrive at the equation 



G%jT{Tj — 1, 

where M is the dimension of the parameter space and the matrix G is defined as follows 



(3.18) 



Gi 



1 d 2 C(r) 



C(0) dndrj 



(3.19) 



Equation (3.18) defines an M-dimensional hyperellipsoid which we take as an approximation to the characteristic 
correlation hypersurface of our random field and we call the correlation hyperellipsoid. This approximation is helpful 
in establishing upper limit estimates of the number of elem entar y cells in the parameter space. The M-dimcnsional 
Euclidean volume V ce n of the hyperellipsoid defined by Eq. ( |3. 18 ) equals 



V ce ii — 



r M/2 



r(M/2 + l)VdetG 



(3.20) 



where L denotes the Gamma function. 

We estimate the number N c of elementary cells by dividing the total Euclidean volume Vtotal of the parameter space 
by the volume V^ G n of the correlation hyperellipsoid, i.e. we have 



Vtotal 
Kell 



(3.21) 



We approximate the probability distribution of J-(£) in eac h cell by the probability distribution po(J-) when the 
parameters are known [in our case it is given by Eq. (3.12a)]. The values of the statistics T in each cell can be 
considered as independent random variables. The prob ability that T docs not exceed the threshold T Q in a given 
cell is 1 — Pf(^F ), where Pp(J r ) is given by Eq. (3.13). Consequently the probability that F does not exceed the 
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threshold T a in all the N c cells is [1 - P F (J r )] Afc 
given by 



The probability Pj, that T exceeds T Q in one or more cell is thus 



Pj(F ) = l-[l-P F (F )} 



(3.22) 



This is the false alarm probability when the phase parameters are unknown. The expected number of false alarms 
Np is given by 



N F = N c P F {To). 
By means of Eqs. ( 3.13 ) and ( |3.21 ), Eq. ( 3.23| ) can be rewritten as 

Vtotal , 



V C eU 



(1 + T ) exp(-^ ). 



(3.23) 



(3.24) 



Using Eq. ( 3.23 ) we can express the false alarm probability Pj, from Eq. ( 3.22 ) in terms of the expected number of 
false alarms. Using lim T1 _ foo (l + = exp(x) we have that for large number N c ^S> 1 of cells 



Pj{To) « 1 - cxp(-7V F ). 



(3.25) 



When the expected number Np of false alarms is small, Np <C 1, we have Pp{T ) « A^. 

It will very often be the case that the filter we use to extract the signal from the noise is not optimal. This may be 
the case when we do not know the exact form of the signal (this is almost always the case in practice) or we choose 
a suboptimal filter to reduce the computational cost and simplify the analysis. For such a case in Paper III we have 
developed a theory of suboptimal filtering. The suboptimal statistics -F S ub has the same form as the optimal statistics 
T and the data analysis procedure consists of maximizing the value of the suboptimal statistics. We have that 



Ei{J-" su b} — 1 



1 



ib) 



(3.26) 



where Ej is the expectation value when the signal is present and the parameter d su ^ is the suboptimal signal-to-noise 
ratio. 

The formulas for the false alarm and detection probabilities have the same form as in the case of the optimal filter 
except that the optimal signal-to-noise ratio is replaced by the suboptimal one. Also for the suboptimal filter the 
number of cells N sc may be different than for the optimal one. Thus the false alarm probability is given by 



Pj F (T ) = l-[l-Pp{T o )] 
The expected number of false alarms N s p reads 

N s p=N sc P F (F ). 

The detection probability takes the form 



N 30 



/>oo 

Poidsnb^o) '•= J pi{d suh ,T)dT, 



(3.27) 
(3.28) 

(3.29) 



where the probability density functional is given by Eq. ( 3.12b| ). Equation ( 3.26| ) implies that maximizing expectation 
value of the suboptimal statistics when the signal is present is equivalent to maximizing suboptimal signal-to-noise 
ratio. We can identify the maximum d smax of d su b as the signal-to-noise ratio of the suboptimal filter. It was shown 
in Paper III that 



= FFrf, 



(3.30) 



where FF is the fitting factor introduced by Apostolatos |n]]. The fitting factor FF between a signal h{t;6) and a 
filter h'(t; 9') (6 and 0' are the parameters of the signal and the filter, respectively) is defined as 



FF 



max 
e' 



(h(t;O)\h'(t;0')) 



(h{t-e)\h{t-e))J(h'{t-e')\h'{t-e')) 



(3.31) 



The fitting factor is a good measure of the quality of the suboptimal filter however th e per forma nce o f the filter can 
only be fully determined by the false alarm and detection probabilities given by Eqs. (3.27) and (3.29). 
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IV. THE FREQUENCY DOMAIN DATABASE 



The data collected by the detector constitute a time series with a sampling interval At. This time domain sequence 
can also be stored in the frequency domain without any loss of information. Organization of the data in the fre- 
quency domain in a suitable way provides a very flexible database that is very useful both for search of the data for 
gravitational-wave signals and for characterization of the noise of the detector. The first step to construct frequency 
domain database consists simply of taking FFT of 2N data points and storing the first N points of the FFT. To 
improve the quality of the time-domain data that will eventually be extracted from the database the original data 
are windowed and overlapped. Each FFT is then normalized and calibrated where the calibration process takes into 
account the transfer function of the detector so that the units of the FFT are strain/ yTlz. Moreover the quality of the 
data represented by FFT is characterized so that it is possible to choose thresholds for vetoing the data or criteria to 
weight them. Both the calibration and the data characterization information are stored in the header that is attached 
to each FFT. It is very important to choose a suitable number N of data in each FFT. This depends on the type of 
signal search that one wants to perform. In the case of the search for continuous sources which is a subject of this 
work we choose N in such a way the maximum expected Doppler shift due to the motion of the detector with respect 
to SSB is less than the width of one bin. In the case of the EXPLORER data this leads to N = 2 16 . For sampling 
time At = 0.18176 s this means length of data for each FFT of around 2/3 of an hour. Another important criterion 
for the choice of the length of data interval for each FFT is that within the interval the noise is stationary. In the 
case of the EXPLORER detector stationarity is in general preserved over 2/3 of an hour intervals. 

The detailed steps to create the frequency domain database are as follows. 

• Each FFT is computed using 2N data, sampled with sampling time At. The data are windowed, in the time 
domain, before the Fourier transform. We use a Hamming window, that is the data yi are multiplied by 
Wi ~ A — B cos i + C cos 2i, where i = {0,2N - 1) • 2n/{2N~ 1), A = 0.54, B = 0.46, and C = 0. 

• The FFTs are calibrated, normalized, and stored in units of strain/ VHz so that their squared modulus is the 
spectrum. 

• The basic FFTs of the database overlap for half their length. The time duration of each FFT is to = 2NAt, 
and a new FFT is done after time to/2. This is important since it avoids distortions in the final time domain 
sequence — this is the well-known "overlap-add" method, described in many data analysis textbooks. 

Once we have a database we can extract from it a time domain sequence of time duration 2MN At and bandwidth 
Av by the following procedure. 

• Take the data from n! = NAv/B bins in the frequency band Av of the actual search, where B is the total 
bandwidth of the detector. 

• Build a complex vector that has the following structure: 

— the first datum equal to zero; 

— the next n' data equal to those from the selected bins of the FFT; 

— zeroes from bins n' + 1 up to the nearest subsequent bin numbered with a power of two; let us say that in 
this way one has obtained n bins; 

— zeroes in the next n bins; one thus ends up with a vector that is 2n long. 

• After the bandwidth has been selected, the data (still in the frequency domain) should be windowed, to avoid 
edge effects in the transformed data. 

• Take the inverse FFT of the vector obtained above. This is a complex time series that is called the analytic 
signal because by construction above its spectrum is zero for negative frequencies. The signal has the bandwidth 
Av that is shifted towards the lower frequencies and the signal is sampled at a sampling rate lower by a factor 
N / n compared to the original time data.[] The time of the first sample here is exactly the same as the time of 
the first datum used for the database and the total duration is also that of the original time stretch. There are 
fewer data because the sampling time of the original data is shorter. 



1 The construction of the analytic signal is a standard procedure of lowpass filtering for a bandpass process. By the fact that 
the analytic signal is zero on the left frequency plan one avoids aliasing effects in the lowpass sampling operation; see, e.g., 
Ref. @. 
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• Remove the window w used when constructing the FFT database simply by dividing the new time domain data 
by w. This operation recovers the original (subsampled) time data because the only regions where the division 
might not work are the edges of the data stream where the value of the function w may be zero, depending on 
the used window (a problem which, of course, has been overcome by overlapping of the FFTs). 

• Repeat the steps outlined above for all the M FFTs. 

• If any FFT under consideration is vetoed or if it is missing the corresponding data are set to zero. 

• Each new group of time domain data is appended to the previous groups after elimination of the overlapped 
data. Since the overlapping involved half the data we eliminate 1/4 of the data at the beginning and end of 
each stream. The first 1/4 data in the first FFT and the last 1/4 in the last FFT can be discarded. The data of 
missing or vetoed periods, set to zero as explained above, are appended to the other stretches in the same way. 

Thus by the above procedure we have a subsampled time domain data stream which represents the analytic signal 
associated with the original data. 

V. POSITION AND VELOCITY OF THE DETECTOR WITH RESPECT TO THE SOLAR SYSTEM 

BARYCENTER 

A. Dealing with time scales 

The data acquisition system of the detector is synchronized to the Coordinated Universal Time (UTC) as dissemi- 
nated by international time services. Thus it is assumed that each datum corresponds to a given UTC. The UTC scale 
is essentially uniform, except for occasional 1 s steps introduced internationally to compensate for the variable Earth 
rotation. When these steps are taken into account, the UTC is reduced to the International Atomic Time (TAI) scale, 
which is uniform and differs from the Terrestrial Dynamical Time (TDT or TT, normally used to describe celestial 
phenomena by astronomers) only by a constant term. In principle, the time argument of positions of celestial objects 
is the Barycentric Dynamical Time (TDB), however it differs from the TDT only by a very small additive term (less 
than 0.001 s in magnitude), thus for all practical purposes they do not have to be distinguished. So we have: 

TDB « TDT = TAI + 32.184s, 

where TAI is obtained from UTC by removing the time steps.^] One can thus relate given UTC with barycentric 
positions of all the major celestial bodies of the solar system. 

However, to relate the position of a point on the Earth to the barycenter of the Earth (the latter being obtained 
from the solar system ephemeris) one has to use yet another time scale — the rotational time scale UT1, which is 
nonuniform and is determined from astronomical observations. The difference UT1 — UTC, which is maintained 
within ±0.90 s, is taken from the IERS tabulations of daily values^] 

B. Topocentric coordinates 

To be able to relate a point on the Earth surface to the solar system barycenter it is necessary to know orientation of 
the Earth in space. The primary effects that should be taken into account are: diurnal (variable) rotation, precession 
and nutation of the Earth rotational axis, and polar motion. The precession and nutation can be accounted for 
by applying standard astronomical theories. The remaining two effects are unpredictable for a longer future, so 
observational data must be used.^ 

The polar motion is taken into account by modifying the conventional geographical coordinates of a point on the 
Earth: 



2 The time steps are available in the form of a table; see http://hpiers.obspm.fr/webiers/general/earthor/utc/ 
tablet .html. 

3 They are available as eopc04.yy files, where yy stands for a two digit year number (e.g. 99 for 1999, and 00 for 2000); see 
http : / /hpiers . obspm . f r/ iers/ eop/eopc04. 
4 For past years, since 1962, the data necessary for reduction are included in eopc04.yy files mentioned in footnote til 
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<fi = <fio + P x cos A — Py sin A , (5.1a) 

A = A + (P x sin A + P v cos A ) tan<^ , (5.1b) 

where </> is the conventional geographical latitude, A — the conventional longitude, and P x and P y are the coordinates 
of the pole with respect to the Conventional International Origin. 

The rotational angle of the Earth is included through conversion of UTC to UT1 (the quantity UT1 — UTC is 
tabulate d in the IERS files). UT1 serves to calculate the apparent sidereal time 9, which is essentially equal to cf> r + U r t 
of Eqs. fl2.3| ) and ( |2.5| ) plus the nutation in longitude projected on the celestial equator. The sidereal time in turn is 
used to find rectangular coordinates of the point (the detector) in the IERS celestial reference frame: 

xe = rcosO, (5.2a) 

He = rsin#, (5.2b) 

ze = 6sin-0 + /isin0, (5.2c) 

where 

r = a cos ip + h cos 4> (5-3) 

is the equatorial component of the radius vector, tp — arctan(6tan <j)/a) is the reduced latitude, h is the height above 
the Earth ellipsoid, and a — 6378.140 km and b = a(l — 1//) (where / = 298.257) are the semiaxes of the Earth 
ellipsoid. These coordinates are affected by the polar motion through dependence of r and ip on (f>, and 9 on A. 

The polar motion (P x and P y ) and UT1 — UTC quantities are linearly interpolated between the daily IERS values. 

Since these coordinates are naturally referred to the epoch of date, they are further precessed back to the standard 
epoch J2000. The precessed Cartesian coordinates (xe,Ue, 2e)20oo may now be straightforwardly added to coordinates 
of the Earth barycenter with respect to the solar system barycenter (which are described in the next paragraph) to 
obtain the desired position vector of the detector. 

C. Barycentric coordinates of the Earth (JPL ephemeris) 

For computing the coordinates of the Earth barycenter, relative to the SSB, use is made of the latest JPL Plane tary 
and Lunar Ephemerides, "DE405/LE405" or just "DE405"J3 created in 1997 and described in detail by Standish [p_3| - 
It represents an improvement over its predecessor, DE403. DE405 is based upon the International Celestial Reference 
Frame (ICRF; an earlier popular ephemeris DE200, which has been the basis of the Astronomical Almanac since 
1984, is within 0.01 arcseconds of the frame of the ICRF). It constitutes of a set of Chebyshev polynomials fit with 
full precision to a numerical integration over 1600 AD to 2200 AD. Over this interval the interpolating accuracy 
is not worse than 25 meters for any planet and not worse than 1 meter for the Moon. The JPL package allows a 
professional user to obtain the rectangular coordinates of the Sun, Moon, and nine major planets anywhere between 
JED 2305424.50 (1599 DEC 09) to JED 2525008.50 (2201 FEB 20). 

In the application described in this paper we have used only a 21-year (1990 to 2010) subset of the original ephemeris. 
The ephemeris gives separately the position of the Earth-Moon barycenter and the Moon's position relative to this 
barycenter. The Earth position vector is thus calculated as a fraction (involving the masses of the two bodies) of the 
Moon's one and opposite to the latter. 

Finally, the vector traveled by the Sun towards its apex (with the speed of 20 km/s) between J2000 and the epoch of 

observation is added to the Earth barycentric position. The direction of solar apex is assumed at 18^ in right ascension 
and 30° in declination at the epoch J1900. This direction is precessed to J2000. What is commonly known as the 
solar apex refers rather to solar system barycenter apex. However since the apex motion is known only approximately, 
really there is no need to distinguish between motions of the Sun and of the barycenter. 



5 It is available via the Internet (anonymous ftp: navigator.jpl.nasa.gov, the directory: ephem/export) or on CD (from the 
publisher: Willmann-Bell, Inc.; http://www.willbell.com/software/jpl.htm). 
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D. Velocities 



Although the primary concern in this project is to convert positions, the velocity of the gravitational-wave detector 
relative to the SSB may prove to be useful in future analyzes of spectra obtained from the acquired data. 

The velocity of the detector is the sum of diurnal rotational motion of the Earth, Earth motion in space relative to 
the SSB, and motion of the solar system itself towards the apex. The last of the named components, towards the apex, 
has been described in the previous paragraph. The Earth barycentric velocity vector is obtained directly from the JPL 
ephemeris along with the position vector. Finally, the detector motion relative to the Earth barycenter is represented 
by a vector of constant length v := 2irr/ (sidereal day) directed always towards the east in the topocentric reference 
frame. Thus this diurnal velocity vector has the following Cartesian components (in the barycentric reference frame): 

V x = v o cos(0 + 7t/2), (5.4a) 

V y = w sin(6» + 7r/2), (5.4b) 

V z = 0. (5.4c) 

This vector is rotated back to J2000 by the precessional angle. 

A short description of the FORTRAN code that generates both the position and the velocity of a detector with 
respect to the SSB is given in Appendix |c]. 



VI. A LINEAR FILTER 



The phase of the gravitational- wave signal contains terms arising from the motion of the detector with respect to the 
SSB. These terms consist of two contributions, one which comes from the motion of the Earth barycenter with respect 
to the SSB, and the other which is due to the diurnal motion of the detector with respect to the Earth barycenter. 
The first contribution has a period of one year and thus for shorter observation times can be well approximated by a 
few terms of the Taylor expansion. The second term has a period of 1 sidereal day and to a very good accuracy can 
be approximated by a circular motion. We thus propose the following approximate simple model of the phase of the 
gravitational- wave signal: 

s 

*»(*) = p + p t + ^2p k t k+1 + Acos(tt r t) + B sm(tt r t), (6.1) 
fe=l 

where Q r is the rotational angular velocity of the Earth. The parameters A and B can be related to the right ascension 
a and the declination S of the gravitational-wave source through the equations [cf. Eq. (18) in Ref. 0] 

A = — — cos5cos(a: — <f) r ), (6.2a) 
c 

B = — — cos<5sin(a — <f> r ), (6.2b) 
c 



where ujq is the angu lar f requency of the gravitational- wave signal and r is defined in Eq. (5.3) 



The phase model (6.1) has the property that it is a linear function of the parameters. We have shown in Paper 
III that for linear phase models the optimal statistics is a homogeneous random field and consequently the statistical 
theory of signal detection des cribed in Sec. Ill of the present paper applies to this case. The polynomial in time part 
of the phase [cf. Eq. (|6.l])1 contains two contributions. The first one comes from the intrinsic frequency drift of 
the gravitational waves emitted by a source. For example, if the source is a spinning neutron star, the frequency of 
the gravitational waves it emits can evolve as the frequency of the revolution of the star. In general the star will lose 
its energy and will spin down. This evolution of the frequency can be approximated by a Taylor series. The second 
contribution comes from the Taylor expansion of the motion of the Earth around the Sun. It is clear that the longer 
the observation time of the signal the more terms of the Taylor expansion we need to include in order to accurately 
approximate the true signal. 



In order to verify the accuracy of the linear model (6.1) we have calculated the fitting factor FF between the linear 
model and the accurate model of the signal. As the accurate model of the phase of the signal we have taken the 
following model: 
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8+1 

fc=i 



fc+i 



s+1 



U> 



fc=l 



n • r SSB (t) 



(6.3) 



where is the signal's frequency too shifted by a known frequency u> s , i.e. ojd = loq — lo s . The down-conversion 
frequency lu s (as well as the bandwidth of the signal) can be chosen arbitr arily and an appropriate time sequence 
can be extracted from the frequency domain database as described in Sec. |y[. The parameters lju are spin-down 
parameters and they arise by approximating the intrinsic frequency evolution of the sign al by a Taylor expansion. In 
the accurate model (^^) we include one more spin down than in the linear model ( |6.l| ). This additional spin down 
serves to represent the uncertainty of our model. 



For the case when the signal is narrow-band around some frequency cvq the formula (3.31) for the fitting factor can 
be approximated by 



FF a max (cos [* a (i; £J - *,(t; fl]) , 



where C a = (^o, <*>i, . . . , w s +i, a, 5) a re t he parameters of the accur ate m odel (S.3) and £ = (p,po,Pi, 
are the parameters of the linear model (6.1). The linear phase model (6.1) can shortly be written as 



s + 3 



(6.4) 

,Ps,A,B) 

(6.5) 



i=0 



The phase in the accurate model ( p.3j ) can be expressed as a sum similar to the sum from Eq. ( |6.5| ) plus a certain 
remainder r that we assume to be small: 



s+3 



*-(*;Ca) = £C»(C,) *<(*) + K*;Ca)- 



(6.6) 



i=0 



Note that the coefficients £ ai are some functions of the accur ate phase parameters £ a . In numerical calculations 
described below the vector tssb in the accurate phase model ( |6.3| ) was computed by our Top2Bary code described 
in Appendix C, and then the coefficients £ ai and the residual r(t; C a ) were computed by making a least-squares fit 
(within the observatio nal i nterval) of the templa te T2 i (, ai h(t) to the accurate phase ^ a (t; C a )- 
Making use of Eqs. (6J5) and (3.6), from Eq. (5.4 ) one gets 



FF: 



max ( cos 
C \ 



s+3 



J2(C ai -Ci) h(t)+r(t) 



max 1 
C 



s+3 



t,J=0 



s+3 



cosr(t) - ^ (C ai - Ci) sinr(t) ) , 



(6.7) 



i=0 



where the last approximation was obtained by Taylor expansion of the cosine function and by keeping only the 
quadratic terms in the (assumed to be small) di fferen ces C, ai — Q. It is now easy to find in Eq. (6.7) the maximum over 
the parameters £ as the right-hand side of Eq. (6/7) is quadratic in these parameters. The values of the parameters 
that maximize the fitting factor are given by the solution of the following set of s + 4 linear equations: 



s+3 



where 



EMc QJ -G) = - L > 

3=0 



Ltj := (k(t) I j (t) cos r(t)) 



Li 



(h(t) smr(t)} . 



(6.8) 



(6.9a) 
(6.9b) 



We have used t he ab ove prescription t o ca lculate the fitting factor numerically. Once we have found the parameters 
of the extremum ( |6.7| ) by solving Eqs. ( |6.8[ ), we have used these values as input to an optimization routine (b ased 
on Nelder-Mead simplex algorithm) to find an accurate value of the fitting factor directly from the formula (6.4) 
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TABLE I. Adequacy of the monochromatic linear model. 
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without the Taylor expansion. In our calculation we have used the following estimates of the maximum values of the 
spin-down parameters 

Klmax = I , (6.10) 

where r m i n is the minimum characteristic spin-down age of the neutron star. These estimates were adopted in the 
previous papers of this series and they were taken from the work of Brady et al. jlj . We have computed the fitting 
factor for three phase models with s — 0,1,2, i.e. for a monochromatic signal, 1-spindown signal, and 2-spindown 
signal. We have carried out computations for two values of the spin-down age (r m j n = 40 years and r m j n = 1000 years), 
and for different values of the signal's frequency /o = wn/(27r) and the observation time T Q . For each case we have 
made calculation for a grid of positions of the source in the sky. We have chosen our observations to start from the 
beginning of the year 2000 (Julian day = 2451545.0) and the position of the detector to coincide with the geographical 
location of the EXPLORER resonant bar. 

Our results are summarized in Tables |, |l[ and III . For each grid of positions of the source in the sky we have found 



the minimum value FF m j n of the fitting factor (the worst case). In the tables we have given the maximum length 
of the observation time for which FF m i n is greater than 0.9, 0.9 1 / 3 , and 0.999. The conservative value of the fitting 
factor equal to 0.9 1 / 3 ~ 0.967 comes from the arguments of Apostolatos jll| by which such a fitting factor leads to 
affordable 10% loss of events. The ultraconservative value of the fitting factor equal to 0.999 gives a negligible loss of 
0.3% of events. 

From the results presented in Tables GL [n], and III it follows that the monochromatic linear model is adequate 



for a few hours of the observation time, 1-spindown model for a few days of the observation time, and 2-spindown 
model for around 1 week of the observation time. For several months of the observation time we do not expect 
to fit an adequate linear model however we know that for such long observation times a coherent all-sky search is 
computationally prohibitive Jl|,||. Thus we conclude that in realistic (i.e. computationally feasible) coherent searches 
of not more than 1 week duration a satisfactory linear model can be chosen. 



VII. SPACING OF FILTERS, THE SEARCH ALGORITHM, AND COMPUTATIONAL REQUIREMENTS 



A. Grid of templates 



In this section we shall present a construction of a grid in the parameter space on which the statistics T will be 
calculated in order to search for signals. We assume that as filter (or template) that we use in order to calculate the 
statistics we shall use an approximations of the signal by a suitable linear model studied in Sec. VI. We would like 
to choose the grid in such a way that there is no loss of potential gravitational-wave signals. In order to determine 
an appropriate grid of templates we shall use the correlation function of the statistics T . We shall choose the grid in 
such a way that the correlation function for two filters evaluated for parameters at two neighboring points of the grid 
is not less than a certain specified value. 

We thus consider here a constant amplitude signal 
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TABLE II. Adequacy of the 1-spindown linear model. 
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TABLE III. Adequacy of the 2-spindown linear model. 
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h (t; ho, * , = sin [* (i; £) + $„] , 



(7.1) 



where $o is the initial phase of the wave form and the vector £ collects all other phase parameters. For the signal 



(7.1) the optimal statistics T can be written as 

2T n 



S h (fo) 



x(t)cos$(t)) 2 + (x(t)sm<f>(t)) 2 . (7.2) 



In Sec. IV of Paper III we have shown that when the noise in the detector is a zero mea n stationary Gaussian process, 
then the correlation function C [defined in Eq. ( 3.1 5| ) ] of the random field J 7 , Eq. fl7.2| ), can be approximated by the 



formula 

C{t €') « <cos[$(<; i) - *(f; £')]) 2 + <sto[3(*5 £) - *(t; £')]> 2 . (7.3) 



When the phase $ is a linear function of the parameters the correlation function C of Eq. ( |7.3| ) depends only on 
the difference r — and can be written as 

C(t) « (cos[$(i; r)]> 2 + (sin[$(t; x)]) 2 . (7.4) 

The fact that for linear phase models the correlation function C depends only on the difference £ — £' between the 
parameters £ of the signal and the parameters £' of the filter, and not on their absolute values, implies that the 
spacing between the filters will be independent of the values of the parameters. This is one of the advantages of the 
linear filter models. 



We further approximate the correlation function C of Eq. (7.4) by Taylor expansion around t = 0. Keeping terms 
at most quadratic in tj one gets 

C(r)«l-^f 4J ^, (7.5) 

where the matrix T has the components 

r y :=/f£fiW|£V|t). (T.6) 



dn drj I \ dn / \ dr 3 



The matrix T is the reduced Fisher information matrix for the signal ( |7.1| ) where the initial phase parameter $ nas 
been reduced (see Appendix B of Paper III for details of this reduction procedure) . 

The statistics T can be expressed as a Fourier transform where the parameter po corresponds to angular frequency 
and consequently T can be calculated using the FFT algorithm. However the FFT gives the values of the statistics on 
a certain grid of frequencies called Fourier frequencies. These frequencies in normalized units are separated by factor 
of 2tt. Thus when the true frequency falls between the two Fourier frequencies we cannot achieve the theoretical 
maximum of the optimal statistics T . The worst case is when the frequency falls half way between the Fourier 
frequencies. One can easily find that in such a case the signal-to-noise ratio calculated approximately by means of the 
Taylor expansion of the correlation function is equal only to 0.42 of the optimal signal-to-noise ratio. This would lead 
to a drastic loss of signals. Therefore we need a finer grid. A way to achieve this and still take advantage of the speed 
of FFT is to pad the time series with zeros.[] Padding with zeros of the length of the original time series gives a grid 
that is twice as fine as the Fourier grid i.e. the difference between the Fourier frequencies is equal to ir. Then in the 
worst case the signal-to-noise ratio is 0.89 of the optimal one. As we shall see later from Monte Carlo simulations this 
choice of spacing of the frequencies leads to a search algorithm where there is no loss of events and rms errors of the 
estimators of parameters arc close to the theoretical minimum. It is also useful to note that the signal-to-noise ratio 
calculated form the exact correlation function formula yields the values of the fractions of the optimal signal-to-nose 



Padding time series with zeros of the length of the original time series in real signal search means that we would need to 
calculate FFTs of twice the length of the original data. This means doubling the computational time and the computer memory 
used. To avoid this pulsar astronomers have invented special interpolation algorithms that work in the Fourier domain and 
give twice as fine Fourier grid from the FFT of the original data. In the analysis of EXPLORER data we shall use one such 
algorithm provided to us by Duncan Lorimer. 



1G 



ratio equal to 0.63 and 0.90 for the difference between the Fourier frequencies equal to 2tt and 7r, respectively. This 
indicates the limits of the validity of the Taylor expansion. 

We introduce a convenient normalization of the spin-down parameters p k of the linear model ( |3.l| ), namely 



p h :=p k T* + \ fc = 0,. 



(7.7) 



This is equivalent to using a time coordinate t := t/T a normalized by the total observation time T Q (then the normalized 
time duration of the signal is always 1). Using definition (7.7) the linear phase model ( |6.lD , after dropping the initial 
phase parameter p, can be written as: 



J- n 



k=l X 



k+1 



A cos(£l r £) + B sm(fi T .i), 



where s is the number of spin downs included. In the quadratic approximation (|7 
equation 



1 



const, 



(7.8) 

and for the phase model (|7.8]), 

(7.9) 



describes the surface of the (s+3)-dimensional correlation hyperellipsoid. 

It is clear that we cannot fill the parameter space completely with hyperellipsoids. But such filling can be done by 
means of some prisms constructed with the aid of the correlation h yper ellipsoid. We shall illustrate our construction 
in the special case of 1-spindown (i.e., s = 1) linear phase model (7.8). In this case the components of the matrix 
r as functions of the observation time are given by [here the order of the phase parameters is as follows: r = 
(Ap^Ap^AAAB)] 



1/12 
1/12 




1/12 
4/45 
l/(2n 2 7r 2 ) 



\ -l/(2mr) -l/(2n7r) 



-l/(2n7r) \ 

l/(2n 2 7r 2 ) -1/(2ti7t) 
1/2 
1/2 / 



(7.10) 



where n is the observation time expressed in sidereal days, i.e., n := (Q r T Q ) / (2it) . 

In the first step we choose spacing in the frequency parameter p to be equal t o n. Then we compute the value Co 
of the correlation at the surface of the 4-dimensional correlation hyperellipsoid (7.9) by r equ esting that the p axis 
intersects this surface at the points Ap — ±n/2. Substituting r = (Ap , 0,0,0) into Eq. ([7j|), one gets 



Co 



7T" 

48 



0.79. 



(7.11) 



Making use of Eqs. ( 7.1C ) and ( |7.1l| ), the general equation (7.9) can be written in the case of 1-spindown signal in 
the form [we also substitute r = (p ,p 1 , A, B)] 



A 2 +B 2 



where we have introduced 



lfo + lpoPi + -^Pl + 2( 2 Api 



2(B(p a 



-*>-24=°' 



1 

777T 



(7.12) 



(7.13) 



We first construct the elementary cell of the grid in the filter space (which is the space spanned by the parameters 
p lt A, and B). We consider the p = cross section of the 4-dimensional hyperellipsoid ( 7.12| ). This cross section 
defines the 3-dimensional ellipsoid in the (p l5 A, B) space: 



(A + C^xf + iB-CPtf = 



as ^2 



24 



45 



(7.14) 



(let us note that ^ — Q 2 — C 4 > for n > 1). Next we take the cross sections of the ellipsoid (7.14) with the planes 
p 1 = ±5 for some positive constant 5. These cross sections define two adjacent circles. We inscribe two regular 
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hexagons into these circles. The hexagons form the bases of the inclined prism inscribed in the ellipsoid (7.14). The 
3-dimensional volume of this prism reads 



V(S)=3y/3 



— - — - c - c 

24 V 45 



6. 



Then we maximize the volume V(S) with respect to 5. The function V(S) attains its maximal value for 

7T 



6V2\/--C 2 -C 4 



(7.15) 



(7.16) 



The maximal value ^(^max) of the volume ( 7.1 5| ) is the volume of the elementary cell in the filter space. It reads 



(7.17) 



2.05, 1.35, 



Equation (7.17) gives the following values of the volume V gv of the elementary cell in the filter space: V gI 
1.29, 1.27 for n = 1, 2, 3, 4, respectively. 

The elementary cell in the 4-dimensional space (p Q ,p 1 , A, B) we construct as follows. The 3-dimensional prism of 
maximal volume, described above, which lies in the p = subspace, we parallelly translate in four dimensions, in 
two opposite directions, into the subspaces p Q — — ir/2 and p Q = +ir/2. As the direction of translation we choose one 
of the principal directions^ of the hyperellipsoid ( 7.12 ). Such constructed elementary cell is the 4-dimensional prism 
which bases are two adjacent 3-dimensional hexagonal prisms lying in the p = — n/2 and p = +ir/2 subs paces . 

It is clear from the construction that our elementary cell will stick out the 4-dimensional hyperellipsoid ( [7.12 ) . For 
the case of n — 2 we have calculated the correlation function C(t) = 1 — J2» j ^ij T i T j a ^ a ^ vertices of the elementary 
cell and we have found that the smallest value of C(r) is equal to 0.77. Thus the grid constructed above ensures that 
the correlation between the filter and the signal is not less than 0.77. The Monte Carlo simulations presented below 
show that using such a grid we do not lose any signals and that rms errors of the parameter estimators are very close 
to the minimum ones allowed by the Cramer-Rao bound. 

If as a base of the elementary cell in the filter space we choose a square instead of a regular hexagon the volume of 
the cell (independently of the value of n) decreases by a factor of 3-^/3/4 ~ 1.3. 



B. Monte Carlo simulations 

In order to test the effectiveness of the chosen grid we have made Monte Carlo simulations of the detection of our 
signal in the noise and estimation of its parameters. In the simulations we have taken the signal s to be 1-spindown 
linear model with a constant amplitude h$: 



s(t) = ho exp < i 



P + Pot+Pit + A cos 



B sin 



/ 2Tint\ 



V T J 



(7.18) 



where T a is the obser vation time and n is the integer equal to the number of sidereal days of observation in real search. 
In the case of signal ( 7.18j ) the maximum likelihood detection involves finding the global maximum of the functional 
T s with respect to the parameters po, pi, A, B. The functional T s is given by 



2|AT 
ShT ' 



(7.19) 



where 



7 The reasonable result one obtains only when translating along this principal direction of the hyperellipsoid (7.12), which 
almost coincides with the p axis. 
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X 



x(t) exp 



Pi t + A cos 



/ 2nnt 



V T J 



Ba "{— ) 



exp(-ipot) dt. 



(7.20) 



In Eq. ( 7.20| ) the data x(t) = s(t)+n(t), where n(t) is the stationary noise with spectral density Sh- In our simulations 
we have generated white and Gaussian noise n(t). Thus our statistics J- s consists of taking the modulus of the Fourier 
transform of the data demodulated for a grid of parameters p\, A, and B. To find the maximum of T s we have 
used a hierarchical algorithm consisting of two steps: a coarse search and a fine search. The coarse search involves 
calculation of T s on the grid in the parameter space constructed in Sec. VII A and finding the maximum value of T s 
on that grid. The values of the parameters of the filter that give the maximum are coarse estimates of the parameters 
of the signal. The fine search involves finding the maximum of T s using a maximization routine with the starting 
values of the parameters equal to the coarse estimates of the parameters. As a maximization routine we have used 
the Nelder-Mead simplex algorithm. This algorithm involves only the calculation of the function to be maximized at 
certain points but not its derivatives. We plan to use the above hierarchical procedure in analysis of real data from 
the EXPLORER detector. 

The procedure described above differs from another two-step hierarchical algorithm proposed by Dhurandhar and 
Mohanty The first step of the two procedures is the same but in the second step Dhurandhar and Mohanty 

propose to use a fine grid in the parameter space around the maximum obtained from the coarse search whereas we 
propose to use a maximization routine. 



We have performed Monte Carlo simulations by generating a signal s(t) given by Eq. (7.18) and adding white 
Gaussian noise n(t). By adjusting the amplitude hg we have generated a signal with a chosen optimal signal-to-noise 
ratio d. In our simulations we have chosen the observation time to be n — 2. We have done 1000 runs for several 
values of d. We have compared the standard deviations of the parameter estimators obtained from the simulations 
with the theoretical ones calculated from the Cramer-Rao bound. We have also compare d the probability of detection 
of the signal obtained from the simulations with t he t heoretical one calculated from Eq. (3.14). In our simulations we 
have used the hexagonal grid constructed in Sec. |VIl| A. 

The results of our computations are depicted in Fig. |]. We have observed that above a certain signal-to-noise ratio 
(that we shall call the threshold signal-to-noise ratio) the results of the Monte Carlo simulations agree very well with 
the calculations of the rms errors from the covariance matrix. However below the threshold signal-to-noise ratio they 
differ by a large factor. This threshold effect is well-known in signal processing [ fl6[ and has also been observed in 
numerical simulations for the case of a coalescing binary chirp signal JlJJl^]. As was explained in Sec. VII of Paper III 
this effect arises because sometimes the global maximum of the functional T s occurs as a result of noise and not the 
signal. This happens the more often the lower the signal-to-noise ratio. Following the theory of this effect developed 
in Paper III we have calculated the approximate rms errors of the estimators of the parameters. They are shown as 
thick lines in Fig. [lj 



The comparison of probability of detection obtained from the simulations with the theoretical formula (3.14) shows 
that for lower signal-to-ratios we have more detections than expected. This is because the theoretical formula assumes 
that when the signal is detected its parameters are located in the cell corresponding to the true parameters of the 
signal. However in practice for lower values of d as a result of the noise the signal may drift to neighboring cells. 
Above the threshold signal-to-noise ratio we find that simulated probability of detection almost exactly agrees with 
the theoretical one and that we are not losing any signals. 



C. Estimation of parameters of the signal 



Once the optimal statistics, calculated with the linear filter of Sec. VI, crosses a chosen threshold we register the 
estimates of the parameters po ; Px-, A, and B. We would like to estimate physically interesting parameters: frequency, 
spin down , declination and right ascension of the source. We propose the following algorithm to achieve this. From 



Eqs. (7.21) and least-squares fit of the linear model to the accurate model of the signal we calculate the approximate 
values of parameters u>q, u)\, S, a. Then using the accurate model of the signal as a filter we search for accurate 
estimates of the parameters on a small grid around the approximate estimates. 

We have made Monte Carlo simulations of the above procedure by making 100 runs for one position of the source 
in the sky and for five signal-to-noise ratios: 8, 12, 16, 20, 24. First of all we have found that probability of detection 
of the signal agrees very well with the theoretical one and thus we are not losing any signals. Concerning accuracy 
of the parameter estimation we have obtained that rms errors for declination and right ascension were very close to 
their Cramer-Rao bounds, however the rms errors for frequency and spin down were worse. The mean biases in the 
estimators of the parameters ujq, u>\, S, and a were less than 0.1%, 0.5%, 0.005%, and 0.01%, respectively. 
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Frequency 



Spin down 






FIG. 1. Monte Carlo simulations of the rms errors of the parameter estimators of the signal given by Eq. ( 7.18 ) for n = 2. 
The x axes are labeled by the optimal signal-to-noise ratio d. The y axes are labeled by the standard deviation a except for 
the case of the amplitude parameter where the relative rms error oy := a /ho is given. The results of simulations of 1000 runs 
for each value of d are marked by the circles. The thin solid lines are calculated from covariance matrix and they constitute 
the Cramer-Rao bound, i.e. the smallest rms error for an unbiased estimator. The thick solid lines are calculated from a model 
that takes into account the false alarms. 
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D. Amplitude modulation 



Approximation of the continu ous gravitational- wave signal by a simple linear model and consequently approximation 



of the optimal statistics, Eq. (3.9), by a homogeneous random field were studied under the assumption that the 
amplitude of the signal was constant. However as the modulation of the amplitude of the signal is very slow (it 
changes on a time scale of one sidereal day) in comparison to the rapid modulation of the phase we expec t tha t our 



linear model of the phase is a satisfactory approximation of the phase of the optimal filters given by Eqs. fTlc| ). We 



also expect that the calculation of the number of elementary cells [by means of Eq. ( 3.2l| ); except for a factor of 2, 



see below] and construction of the grid in the parameter space given in Sec. VII A above are valid for the amplitude 
modulated signal. 

To calculate the optimal statistics, Eq. ( |3.9| ), we first need to correct the time series for amplitude modulation i.e. 
multiply the time series by the functions a(t) and b(t). Our grid is on the space parameterized by po, pi, A, and B 
whereas the functions a(t) and b{t) depend on the declination 5 and the right ascension a. From Eqs. (|6j) we have 
the following expressions for 5 and a in terms of A and B: 



V A 2 + B 2 

5 = ± arccos | jj^p ] , (7.21a) 



c 

a = (j) r + arctan (^~^J ■ (7.21b) 
Since for each set of parameters A and B there are two sets of parameters 5 and a the number of cells that we obtain 



from Eq. (3.21) needs to be multiplied by a factor of 2 and this number is to be used to compute from Eq. (3.22) the 
threshold corresponding to a given false alarm probability. 

We see that the equation for declination <5 of the source depends on the angular frequency ujq of the signal that we 
do not know before we detected the signal. The uncertainty in 5 affects the constant amplitudes in front of the time 
dependent modulation factors. However for the case of the EXPLORER detector data (which we plan to analyze 
emp loyin g the methods developed above) we have that frequency /o is around 1 kHz within the band B ~ 1 Hz (see 
Sec. VH]| below) . Thus by choosing in the equation for declination 5 the frequency equal to the middle frequency of 



the band to be analyzed instead of the unknown frequency ojq, we shall lose only a fraction ~ B/(2Jq) ~ 0.0005 of 
the signal-to-noise ratio. 



E. Computational requirements 

To estimate the computational requirements to do the search we calculate the number of FFTs that one needs to 
compute. This number depends on the total volume Vf of the filter space, which is given by 

f f / -^ ma xT 2 /(2r min ) 

V F = dAdB dp^ (7.22) 

J JB 2 {0,r o ) ^-w max T 2 /(2T min ) 

where a-> max is the maximum angular frequency of the signal and #2(0, r„) is a 2-dimensional disc in the (A, B) plane 
of radius r , 

r = (7.23) 

c 

Thus we get 

3 2 T 2 

V F = 9_ ( 7 24 ) 

C Tmin 

The number of grid points on which the optimal statistics J 7 , Eq. ( [3.9] ) , should be calculated is obtained by dividing 
the volume Vp of the filter space by the volume V gI of the elementary cell of the chosen grid in the filter space. In the 
previous subsection we have found that for each set of parameters A and B there are two sets of parameters S and a 



[ see E qs. (7.21)]. Moreover the optimal statistics involves calculations of two FFTs — F a and Fj, [see Eqs. (3.9) and 
( |3.10 )]. Thus the total number TVfft of FFTs is 4 times the number of grid points and it is given by 
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7V FFT=4 ^. (7.25) 

Assuming that the data processing rate should be comparable to the data acquisition rate and that the data 
processing consists only of computing the FFTs the number V of floating point operations per second (flops) can be 
estimated by 

V = 6A^ FFT [log 2 (2Ai/r o ) + 1/2], (7.26) 
where Av = (w max — w s )/(27r) is the bandwidth of the search. 

VIII. A NARROWBAND ALL-SKY SEARCH OF EXPLORER DATA 

The data analysis tools for searching continuous gravitational-wave signals that we have developed in the previous 
sections of the present work can be applied both to the resonant bar and the laser interferometric detectors. We plan 
to apply these tools to the data from the resonant bar detector EXPLORER^ . The detector has already collected 
many years of data with high duty cycle (e.g. in 1991 the duty cycle was 75%). Our primary objective is to carry out 
an all-sky search. It is a unique property of the gravitational-wave detectors that with a single time series one can 
search for signals coming from all sky directions. In the case of other instruments like optical and radio telescopes 
to cover the whole or even part of the sky requires a large amount of expensive telescope time. The directed search 
of the galactic center has already been carried out and limits for the amplitude of the gravitational waves has been 
established @. 

The EXPLORER detector is most sensitive over certain two narrow bandwidths (called minus and plus modes) 
of about 1 Hz wide at two frequencies around 1 kHz. To make the search computationally feasible we propose an 
all-sky search of data of a few days long in the narrow band were the detector has the best sensitivity. By narrowing 
the bandwidth of the search we can shorten the length of the data to be analyzed as we need to sample the data at 
only twice the bandwidth and thus we reduce the computational time. To reduce the parameter space to search we 
restrict ourselves to only one spin-down parameter. We would also like to use the linear filter model as the search 
template. Then as our frequency is around 1 kHz from Table [l| we read that we can consider up to 2 days of coherent 
observation time in order that the fitting factor is greater then 0.9. For the sake of the FFT algorithm it is best to 
keep the length of the data to be a power of 2. Consequently we choose the number of data points to analyze to be 
N = 2 18 . Then for T a = 2 days of observation time the bandwidth Av of the data will be Av = N/(2T ) ~ 0.76 Hz. 
We also choose to analyze the data for the plus mode which has frequen cy ar ound 922 Hz. With these parameters we 
can calculate the number of filters that we need to compute. From Eq. ( 7.25| ) and for the hexagonal grid over the sky 



we get that TVfft ~ 3.7 x 10 s . Assuming that the data processi ng rat e should be comparable to the data acquisition 



rate the computing power required [calculated by means of Eq. ( 7.26 )] is around 7.7 Gflops. If we allow a month for 
off-line proce ssing of data with the above parameters we only need around 250 Mflops of computer power. 

Using E q. (|3.21 ) we can calculate the number of cells N c in the parameter space. Then from the number of cells, 



using Eq. ( [3.27 ) we can calculate the threshold that we need to set. For example choosing the false alarm probability 



to be 1% wc find by inverting Eq. (3.27) and using Eq. (3.26) that for the parameters that we have chosen above the 



threshold signal-to-noise ratio is 8.3. However we know that in the coarse search we are using a suboptimal filter and 
we are losing signal-to-noise ratio. For the parameters of our search we find that in the worst case the fitting factor 
is 0.94. Moreover due to our discrete grid in the parameter space the signal-to-noise ratio can decrease in the worst 
case by an additional factor of V0.78 ~ 0.88. Thus in the worst case the signal-to-noise ratio of our coarse search can 
be a factor of 0.83 of the optimal one. Hence in order not to lose any signals we lower the signal-to-noise threshold 
by a factor of 0.83. By lower ing the threshold we must make sure that the number of false alarms does not increase 



excessively. From Eq. ( 3.28 ) we find that the expected number of false alarms with the lower threshold is ~ 310. 
This is certainly a manageable number that will insignificantly increase the computational time of the total search. 
We can consider lowering the threshold even further in order to have more candidate events that can be later verified 
using longer stretches of data. For example by lowering the threshold by a factor of 0.8 the expected number of false 
alarms is around 1600 which will still be manageable to verify. 



8 The EXPLORER detector is operated by the ROG collaboration located in Italian Istituto Nazionale di Fisica Nucleare 
(INFN); see http: //www.romal . inf n. it/rog/ explorer/explorer .html. 
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The minimum detectable amplitude ho, i.e. the amplitude for which the signal-to-noise ratio is equal to 1, is given 

by 



h = \ 777" i 



(8.1) 



where Sh is the one-sided spectral density of noise. For 2 days of observation time T and spectral density Sh at the 
plus mode equal to 2 x 10 _42 /Hz, the minimum detectable amplitude is 3.4 x 10~ 24 . In order that we have a detection 
with 99% confidence by the calculation above we need a signal-to-noise ratio of 8.3 and thus the amplitude of the 
signal must be around 2.8 x 10~ 23 . Consequently by estimate given in Eq. (2.7) a continuous gravitational-wave signal 
from a neutron star located at a distance of 1 kpc, spinning with period of 2 ms, and with ellipticity e ~ 10^ 5 will be 
detectable. 
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APPENDIX A: DERIVATION OF THE BAR BEAM-PATTERN FUNCTIONS 



We consider here the response of a bar detector to a weak plane gravitational wave in the long wavelength approx- 
imation. Moreover, we assume that the frequency spectrum of the gravitational wave which hits the bar entirely lies 
within the sensitivity band of the detector. Under these assumptions the dimensionless response function h of the bar 
detector can be computed from the formula (cf. Sec. 9.5.2 in Ref. M) 



h(t) = n • H(t)n , 



(Al) 



where n denotes the unit vector parallel to the symmetry axis of the bar, H is the three-dimensional matrix of the 
spatial metric perturbation produced by the wave in the proper reference frame of the detector, and a dot stands for 
the standard scalar product in the three-dimensional Cartesian space. In the detector's reference frame we introduce 
Cartesian detector coordinates (xd,yd, Zd) with the z& axis along the Earth's radius p oint ing toward zenith, and the 
Xd axis along the bar's axis of symmetry. In these coordinates the vector n from Eq. (Al) has components 



n 



(1,0,0). 



(A2) 



In the wave Cartesian coordinate system (x w , j/ w , z w ) (in which the wave travels in the +z w direction), the three- 
dimensional matrix H of the wave induced spatial metric perturbation has components 



H(t) 



/ h+(t) h x (t) 0\ 

M<) -M*) 

\ / 



(A3) 



where the functions h+ and h y describe two independent wave's polarizations. The matrices H and H are related 
through equation 



H(t) = M{t) H{t) M(tf 



(A4) 



where M is the three-dimensional orthogonal matrix of transfo rma ti on fr om the wave coordinates to the detector 
coordinates, T denotes matrix transposition. Collecting Eqs. ( |Al[ )-( |A4] ) together one can see that the response 
function h is a linear combination of the functions h + and h x : 



h(t) = F+(t)h+(t) + F x (t)h x (t), 



(A5) 
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where F + and F x are beam-pattern functions. 

Because of the diurnal motion of the Earth the beam-patterns F+ and F x are periodic functions of time with 
a period equal to one sidereal day. We want to express F + and F x as functions of the right ascension a and the 
declination 6 of the gravitational-wave source and the polarization angle ip (the angles a, 5, and ip determine the 
orientation of the wave reference frame with respect to the celestial reference frame defined below) . We represent the 
matrix M of Eq. (A4) as 



M = M 3 M 2 Mi, 



(A6) 



where Mi is the matrix of transformation from wave to celestial frame coordinates, Mi is the matrix of transformation 
from celestial coordinates to cardinal coordinates and M3 is the matrix of transformation from cardinal coordinates 
to detector coordinates. In celestial coordinates the z axis coincides with the Earth's rotation axis and points toward 
the North pole, the x and y axes lie in the Earth's equatorial plane, and the x axis points toward the vernal point. 
In cardinal coordinates the (x, y) plane is tangent to the surface of the Earth at detector's location with x axis in the 
North-South direction and y axis in the West-East direction, the z cardinal axis is along the Earth's radius pointing 
toward zenith. Under the above conventions the matrices Mi, M2, and M3 are as follows (matrices Mi and Mi are 
taken from Sec. II A of Ref. [9) 



Mi = 



Mi 



M, = 



/ sin a cos ip — cos a sin 6 sin ip — cos a cos tp — sin a sin 5 sin ip cos S sin ip \ 

— sin a sin ip — cos a sin 8 cos ip cos a sin tp — sin a sin 6 cos ip cos (5 cos tp 

\ — cos a cos 6 — sin a cos S — sin 5 / 

/ sin0cos(0 r + Q r t) sin0sin(0,. + Q r t) — cos(p\ 

— SXD.((pr + &rt) COs((/) r + H r t) 

\ cos (p cos((p r + Q, r t) cos0sin(0 r + Q r t) sin0 / 

/ — sin 7 cos 7 \ 

— cos 7 — sin 7 

V o 01/ 



(A7a) 



(A7b) 



(A7c) 



In Eq. ( A7b) <p is the geodetic latitude of the detector's site, f2 r is the rotational angular velocity of the Earth, and 
the phase <p r defines the position of the Earth in its diurnal motion at t = (th e sum <p r + £l r t essentially coincides 
with the local sidereal time of the detector's site; see Sec. VII B). In Eq. (A7c) 7 determines the orientation of the 
bar's axis of symmetry with respect to local geographical directions, it is measured counter-clockwise from East to 
the bar's axis of symmetry. 

To find the explicit formula for F. 
we arrive at the expressions ( 



and F x we have to combine Eqs. (Al)-(A7). After some algebraic manipulations 
) and (2.5). Equivalent explicit formulas for the functions a and b from Eqs. (2.5) 



can be found in Ref. || where different angles describing the position of the gravitational- wave source in the sky and 
the orientation of the detector on the Earth are used.R 



APPENDIX B: THE TIME AVERAGES (a 2 ) AND (b 2 
1. Laser interferometeric detector 



The time averages (a 2 ) and (fe 2 ) entering Eqs. (3.8), for the observation time T a chosen as an integer number of 



sidereal days, for the laser interferometeric detector take the form (here n is a positive integer): 



9 Our functions a and b from Eqs. (2.5) are identical with the functions S. 
provided the following identification of the variables 77, a, 9, 9 used in 
T) -> -(a - 



■ Q r t), a -> 7r/2 - S, 6 -» tt/2 - (p, * -> 7 - tt/2. 



and 5 X from Eqs. (2.90) and (2.91) of Ref. g, 
with our variables Q r , cp r , a, S, (p, 7 is made: 
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<« 2 > 



T =n2TT/Q, r 



16 



sin 27 



9 cos 4 <P cos 4 5 + ^ sin 2 2<f> sin 2 25 + (3 - cos 20) 2 (3 - cos 25) 2 



<6 2 ) 



+ -COS 2 7 



4 cos 2 </> sin 2 25 + sin 2 0(3- cos 25) 2 



(Bla) 



- 1 - sin 2 27 [(3 - cos 20) 2 sin 2 5 + 4 sin 2 20 cos 2 5 



+ - cos 2 27 (1 + cos 2(j) cos 25) . 



(Bib) 



We see that (a 2 ) and (& 2 ) depend only on one unknown parameter of the signal — the declination 5 of the 
gravitational-wave source. They also depend on the latitude 01 the detector's location and the orientation 7 of the 
detector's arms with respect to local geographical directions. 



2. Resonant bar detector 



For the resonant bar detector the time averages (a 2 ) and (& 2 ), for the observation time T D being an integer number 
of sidereal days, have the form: 



(a 2 ) 



T a =n 27r/O r 



j ( 1 — 3 sin 2 7 cos 2 0) 2 cos 4 5 + i sin 2 7 cos 2 ( 1 — sin 2 7 cos 2 0) sin 2 25 



T =n 2-!r/n r 



+ 8 
1 



sin 2 27 sin 2 + (cos 2 7 — sin 2 7 sin 2 0) (l + sin 2 5) , 
(sin 2 27 cos 2 + sin 4 7 sin 2 20) cos 2 5 



+ i ( cos 4 7 + i sin 2 27 sin 2 + sin 4 7 sin 4 ) sin 2 5. 



(B2a) 



(B2b) 



As in the case of the laser interferometeric detector, the time averages (a 2 ) and (6 2 ) depend on the declination 5 
of the gravitational-wave source, the latitude <p of the detector's location, and the orientation 7 of the detector's axis 
of symmetry with respect to local geographical directions. 



APPENDIX C: THE TOP2BARY PROCEDURE 



The algorithms described in Sec. V were implemented in an easily callable subroutine named Top2Bary consisting 
of about 900 lines of FORTRAN code. The ephemcris files (DE405'90. '10, tai-utc.dat and yearly eopc04.yy) 
required by the program must be kept in the directory \Top2Bary\EphData on the current disk. The procedure 
header is shown below. 

subroutine Top2Bary (C JD_E , clat , clong , height , pve , pvo) 

c Procedure to calculate position (km) and velocity (km/s) of an observatory 
c located at given geographical position (conventional coordinates in degrees 
c and height above the Earth ellipsoid in m) at given Julian day. 
c It requires some ephemeris files, placed in the 'path' directory. 



c Argument description: 

c CJD_E - Julian Day number representing the Coordinated Universal Time 
c If CJD_E is negative it is assumed that it is -1* (Ephemeris JD) . 

c clat - observatory conventional geographical latitude (deg) . 
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c clong - observatory conventional geographical longitude (cleg) . 

c height- observatory height above the Earth ellipsoid (m) . 

c pvo - 6 element array of barycentric vectors (3 for position in km, 

c and 3 for velocity, km/s) of the observatory relative to Earth baryc. 

c pve - 6 element array of barycentric vectors (3 for position in km, 

c and 3 for velocity, km/s) of the Earth barycenter relative to SSB. 

c More important subroutines called: 

c polmot(CJD,Px,Py,UTl_UTC,dpsi,deps) - reads polar motion and UT1 data 

c topobs (DJ1 ,glat ,glong, height , pvo , amst) - returns observer's coord, in 'pvo' 

c earthPV(BJD,pve) - returns Earth barycentric position (J2000 frame) in 'pve' 
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